Sobol G with k = 6

\[G = \prod_{i = 1}^k g_i\]

\[g_i = \frac{|4X_i - 2| + a_i}{1 + a_i}\]

Sensitivity Indices: analytical results

\[V_i = \frac{1}{3(1 + a_i)^2}; \]

\[ S_i = \frac{V_i}{\prod_{i=1}^k \left( V_i + 1 \right) - 1} \]

\[S_{Ti} = V_j\cdot\prod_{i\neq j}^k(1 + V_i)\]

Analytical sensitivity indices calculation

The function for calculating the variances and hence the sensitivity indices are defined below:

options(scipen = 999)

set.seed(100)

# Define the function to calculate V_i
calculate_Vi <- function(a_i) {
  return(1 / (3 * (1 + a_i)^2))
}

# Define the function to calculate S_i
calculate_Si <- function(Vi, V_product) {
  return(Vi / (V_product - 1))
}

# Define the function to calculate V_Ti
calculate_VTi <- function(Vi, V) {
  return(Vi * V)
}

# Define the function to calculate S_Ti
calculate_STi <- function(VTi, V_product) {
  return(VTi / (V_product - 1))
}

# Main script
a <- c(0, 2, 10, 30, 99, 99) # Example values for a
k <- length(a)

# Calculate V_i for each a_i
Vi <- sapply(a, calculate_Vi)

# Calculate the product of (1 + V_i) for total variance
V_product <- prod(1 + Vi)

# Calculate V_Ti for each a_i
VTi <- sapply(1:k, function(i) calculate_VTi(Vi[i], prod(1 + Vi[-i])))

# Calculate first-order sensitivity indices S_i
Si <- Si_analytical <- sapply(Vi, calculate_Si, V_product = V_product)

# Calculate total-order sensitivity indices S_Ti
STi <- STi_analytical <- sapply(1:k, function(i) calculate_STi(VTi[i], V_product))

The arbitrary coefficients are selected. The higher the value, the less significant the input parameter’s impact on the output.

a <- c(0, 2, 10, 30, 99, 99) # Example values for a
# a <- c(6.42,6.42,6.42,6.42,6.42,6.42)
k <- length(a)

The results are then printed:

# Output the results
cat("First-order sensitivity indices (Si):\n")
print(round(Si, 5))

cat("\nTotal-order sensitivity indices (STi):\n")
print(round(STi, 5))

The sum must be equal to 1.

Numerically calculate indices using the sensobol package.

1. Define the sensitivity parameters and the function under inquiry:

set.seed(2)

k <- length(a)
params <- paste("x", 1:k, sep = "")
R <- 10^3
type <- "norm"
conf <- 0.95

sobol_Fun <- function(X, vector) {
  y <- 1
  
  for (j in 1:length(a)) {
    y <- y * (abs(4 * X[, j] - 2) + vector[j])/(1 + vector[j])
  }
  return(y)
}

2. Sensitivity indices from sensobol are computed for each increasing power of 2 with lapply function:

N_values <- 2^(7:15)

results <- lapply(N_values, function(N) {
  mat <- sobol_matrices(N = N, params = params, type = "R")
  y <- sobol_Fun(mat, a)
  ind <- sobol_indices(Y = y, N = N, params = params, boot = TRUE, type = "norm", conf = 0.95, R = R)
  
  list(N = N, 
       Si = ind$results %>% filter(sensitivity == "Si") %>% pull(original), 
       STi = ind$results %>% filter(sensitivity == "Ti") %>% pull(original))
})

3. the results object is converted to two usable data.frame objects, for an easier use:

The convergence plots for each input sensitivity can be plotted:

# Transform the data to be properly plotted
plot_data <- final_results_Si %>%
  pivot_longer(cols = -c(coefficients, Si_analytical),
               names_to = "N",
               values_to = "Si_value",
               names_prefix = "Si_") %>%
  mutate(N = as.numeric(N))

# Create the plot
Si_ggplot <- ggplot(plot_data, aes(x = N, y = Si_value, color = factor(coefficients))) +
                    geom_line() +
                    geom_point() +
                    geom_hline(aes(yintercept = Si_analytical, color = factor(coefficients)), 
                               linetype = "dashed") +
                    scale_x_log10(breaks = unique(plot_data$N)) +
                    labs(title = "Convergence of Si values to analytical values (dashed lines)",
                         x = "Sample size (N)",
                         y = "",
                         color = "Coefficients:") +
                    theme_minimal() +
                    theme(legend.position = "right")
                  
ggplotly(Si_ggplot)

# Transform the data to be properly plotted
plot_data <- final_results_STi %>%
  pivot_longer(cols = -c(coefficients, STi_analytical),
               names_to = "N",
               values_to = "STi_value",
               names_prefix = "STi_") %>%
  mutate(N = as.numeric(N))

# Create the plot
STi_ggplot <- ggplot(plot_data, aes(x = N, y = STi_value, color = factor(coefficients))) +
                     geom_line() +
                     geom_point() +
                     geom_hline(aes(yintercept = STi_analytical, color = factor(coefficients)), 
                                linetype = "dashed") +
                     scale_x_log10(breaks = unique(plot_data$N)) +
                     labs(title = "Convergence of STi values to analytical values (dashed lines)",
                          x = "Sample size (N)",
                          y = "",
                          color = "Coefficients:") +
                     theme_minimal() +
                     theme(legend.position = "right")

ggplotly(STi_ggplot)

Pass the cursor over the dots the read the sensitivity values.

Given the simplicity of the function, the convergence is reached quickly, even with a small sample size. These plots can thus be readily used for more complex functions, where convergence might be slower and more interesting to assess.